A model for the orientational ordering of the plant microtubule cortical array 
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The plant microtubule cortical array is a striking feature of all growing plant cells. It consists of a 
more or less homogeneously distributed array of highly aligned microtubules connected to the inner 
side of the plasma membrane and oriented transversely to the cell growth axis. Here we formulate a 
continuum model to describe the origin of orientational order in such confined arrays of dynamical 
microtubules. The model is based on recent experimental observations that show that a growing 
cortical microtubule can interact through angle dependent collisions with pre-existing microtubules 
that can lead either to co-alignment of the growth, retraction through catastrophe induction or 
crossing over the encountered microtubule. We identify a single control parameter, which is fully 
determined by the nucleation rate and intrinsic dynamics of individual microtubules. We solve the 
model analytically in the stationary isotropic phase, discuss the limits of stability of this isotropic 
phase, and explicitly solve for the ordered stationary states in a simplified version of the model. 
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I. INTRODUCTION 



Most plant cells grow by uniaxial expansion. Estab- 
lishing and maintaining this characteristic anisotropic 
growth mode requires regulatory mechanisms that are 
robust, and, in addition, sensitive to the cell geometry 
A major role in this process is played by microtubules, 
highly dynamic filamentous protein aggregates that form 
one of the components of the cytoskeleton of all eukary- 
otic organisms (see Ref. [1], chapter 16). In growing 
plant cells microtubules are confined to a thin layer of 
cytoplasm just inside the cell plasma membrane. Here 
they form the so-called cortical array, an ordered struc- 
ture formed by highly aligned (bundles of) microtubules 
oriented transversely to the growth direction [2]. This 
structure is unique to plant cells and there is evidence 
that it controls the direction of cell expansion by guid- 
ing the mobile transmembrane protein complexes that 
deposit long cellulose microfibrils in the plant cell wall 
[3]. These cellulose microfibrils are the main structural 
elements of the cell wall, which, for mechanical reasons, 
are also transversely oriented to the cell axis in growing 
cells [4]. An in vivo image of the array and a schematic 
are shown in Figure 1 . In vivo imaging of microtubules la- 
belled with fluorescent proteins in plant cells by several 
groups has shown how the cortical array is established 
both following cell division and after microtubule depoly- 
merizing drug (oryzalin) treatment [2, 3, 5, 6, 7, 8]. In 
these studies microtubules are seen to nucleate at the cor- 
tex and then develop from an initially disorganized state 
into the transverse ordered array over a time period on 
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Figure 1: (a) Fluorescently labelled microtubules in a Tobacco 
BY-2 cell expressing GFP:TUA6. Image courtesy of Jelmer 
Lindeboom, Wageningen University, (b) Schematic represen- 
tation of a plant cell cortical array. 



the order of one hour. The nature of the self-organization 
process by which the specific spatial and orientational 
patterning of this cytoskeletal structure is achieved is as 
yet only partially understood and forms the subject of 
this work. 

An important aspect of the problem is the nature of lo- 
calization of the microtubules to the cortical region. Flu- 
orescence recovery after photo-bleaching (FRAP) experi- 
ments by Shaw et al. [7] showed that the microtubules are 
fixed in space, so any apparent mobility of microtubules 
is due to 'trcadmilling', the process of simultaneous poly- 
merization at one end and dcpolymerization at the other 
end. So, cortical microtubules do not translate or rotate 
as a whole. The same authors also did not detect de- 
tachment or (re) attachment of microtubules to the cell 
cortex, apart from some growing ends of single micro- 
tubules moving out of focus and found no evidence for 
motors working in the cortical array. These experiments 
indicate that the microtubules in the cortical array are 
fixed to the inside of the cell membrane. Electron mi- 
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croscopy has also shown cross-bridges between cortical 
microtubules and the membrane [9] . It is therefore widely 
assumed that there are linker proteins that anchor the 
microtubules through the plasma membrane to the rigid 
cell wall, although their molecular identity is under de- 
bate [10, 11, 12, 13, 14]. Since the cortical microtubules 
are effectively confined to a 2D surface, they can interact 
through 'collisions' that occur when the polymerizing tip 
of a growing microtubule encounters a pre-existing mi- 
crotubule. The resulting dynamical interaction events 
were first characterized by Dixit and Cyr [15] in tobacco 
Bright Yellow-2 (BY-2) cells. They observed three dif- 
ferent possible outcomes: (i) zippering: a growing micro- 
tubule bending towards the direction of the microtubule 
encountered, which occurs only when the angle of inci- 
dence is relatively small (< 40°) (ii) induced catastrophe: 
an initially growing microtubule switching to a shrinking 
state and retracting after the collision, an effect predom- 
inant at larger angles of incidence and (iii) cross-over: a 
growing microtubule 'slipping over' the one encountered 
and continuing to grow in its original direction. 

There are clearly many coupled mechanisms at work 
in this complex biological system contributing to the as- 
sembly and maintenance of this microtubule cortical ar- 
ray structure. We are interested in understanding what 
are the main contributing factors and how their interplay 
leads to the observed oricntational ordering. With this 
aim we develop a coarse-grained model, incorporating all 
the effects discussed above. Our emphasis on the plant- 
specific biological mechanism of the ordering in the cor- 
tical array distinguishes our approach from earlier work. 

Over the years, various models for self-organization 
of cytoskelctal filaments (and polar rods in general) 
have been proposed [16, 17, 18, 19], and the model by 
Zumdieck et al. [20] was applied to the plant cortex. 
However, in each of these models the filaments are as- 
sumed to have rotational and, in most cases, translational 
degrees of freedom. This is inconsistent with the fact that 
the plant cortical microtubules are stably anchored. In- 
spired by the experimental results of Dixit and Cyr [15], 
Baulin et al. [21] were the first to report on a two di- 
mensional dynamical system of treadmilling and colliding 
microtubules. Their focus was on establishing the mini- 
mal interactions needed to generate dynamical alignment. 
Using stochastic simulations they showed that a pausing 
mechanism, whereby a growing microtubule stalls against 
another microtubule until the latter moves away, can in- 
deed lead to ordering. Stalling, however, is not often 
observed in the cortical array. Moreover their model 
lacks dynamic instabilities, i.e. catastrophes, both spon- 
taneous and induced, and rescues, and employs a form 
of deterministic microtubule motion, which is arguably 
unrealistic in view of the observed dynamics. 

The outline of the paper is as follows. In section II we 
formulate our course-grained model starting from a de- 
scription of the dynamics of individual microtubules. We 
then construct the continuity equations that couple the 
densities of growing, shrinking and inactive microtubule 



segments due to the intrinsic and collisional dynamics. 
In the steady state we can reduce the initial set of equa- 
tions to four coupled non- linear integral equations. We 
then perform a dimensional analysis to identify the rele- 
vant control parameter of the system. In section III we 
present the results of our model. We first solve the model 
analytically in the isotropic stationary state. Using a bi- 
furcation analysis we then determine the critical values 
of the control parameter at which the system develops or- 
dered solutions. We interpret these results in terms of the 
physical parameters of microtubule segment length and 
mesh size. Finally, we formulate a minimal model with 
realistic interaction parameters that we can solve numeri- 
cally to obtain all stationary ordered solutions. We close 
by giving arguments for the stability of these solutions. 
The paper concludes with a discussion section. An ap- 
pendix outlines further details of the numerical solution 
technique employed. 



II. MODEL 

A. Description of the microtubules and their 
dynamics 

As described in the introduction we confine the config- 
uration of the microtubules to a 2D plane. Since collision- 
induced zippering events can cause microtubules to bend 
along the direction of preexisting ones, we divide each mi- 
crotubule into distinct segments with a fixed orientation. 
We treat these segments as straight rigid rods. This is 
justifiable since the persistence length l p of microtubules 
is long (~ mm) compared to the average length of a mi- 
crotubule (~ 10/im) and, as mentioned above, adhesion 
to the plasma membrane further inhibits thermal motion. 

Microtubules arc known to be dynamic in that they are 
continually growing or shrinking by (de)polymerization. 
We use the standard two-state dynamic instability model 
of Dogterom and Leibler [22] which assumes that each 
microtubule has a 'plus' end, located on the final segment 
of each microtubule, that is cither growing (labelled by +) 
with speed v + or shrinking (labelled by — ) with speed v~ . 
This plus end can switch stochastically from growing to 
shrinking (a so-called 'catastrophe') with rate r c , or from 
shrinking to growing (a so-called 'rescue') with rate r r in 
a process known as dynamic instability. 

We model the creation of new microtubules with a 
constant, homogeneous, isotropic nucleation rate r n in 
the plane of the 2D model. In vivo nucleation appears 
to occur at the cortex and has been observed to occur 
in random orientations unattached to pre-existing micro- 
tubules [7]. Although microtubules have also been ob- 
served to nucleate from by 7-tubulin complexes binding 
to pre-existing microtubules [2, 23, 24] we ignore this 
possibility for simplicity's sake. By the same token we 
disregard the possibility of the shrinking of microtubules 
at their less active 'minus' end, leading to motion through 
the 'treadmilling' mechanism [25]. The initial segment of 
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Figure 2: Schematic representation of the model interaction. 
The microtubule of interest is drawn in black and other micro- 
tubules that it encounters are in grey. The active segments of 
the black microtubule have an arrow head indicating growth 
or shrinkage whereas inactive segments end in the junction 
with the following segment depicted by a dot. See also the 
description of the parameters in table I. 



each microtubule therefore remains attached to the nu- 
cleation point in our model. 

We call the final segment of a microtubule, which con- 
tains the growing or shrinking tip, active and all the 
remaining ones, which do not change their length, in- 
active (labelled by 0). A cartoon of an an individual 
microtubule according to these definitions is depicted in 
figure 2. When a microtubule collides with another micro- 
tubule and experiences a zippcring event, its active seg- 
ment is converted into an inactive segment, and a new ac- 
tive segment is created alongside the encountered micro- 
tubule. The inverse can also occur: if the active segment 
shrinks to zero length, a previously inactive segment in 
another direction can be reactivated. An induced catas- 
trophe event simply causes the growing active segment 
to become a shrinking one, as is the case for spontaneous 
catastrophes. Finally, a crossover results in the growing 



active segment continuing to grow unperturbed. 

In figure 3 we present the relative probabilities for zip- 
pcring, induced catastrophes and crossovers as a result of 
collisions between microtubules, based on the data pro- 
vided by Dixit and Cyr [15]. We assume that there are 
no microtubule polarity effects, as they were not reported. 
The probabilities P z (0 - 0'), P x (0 - 0') and P c (9 - 0') for 
zippcring, cross-overs and induced catastrophes respec- 
tively are therefore even functions of the angle difference 
6 — 9' defined by their values on the interval [0, ^]. In 
this article we will use only the following minimal set of 
properties, which are qualitatively supported by the data. 
Firstly, zippcring becomes less likely for increasing angle 
of incidence, and is effectively zero at — 6' = J, which 
is reasonable as the energy associated with bending the 
microtubule increases with angle. Secondly, the proba- 
bility for induced catastrophes monotonically increases 
with increasing angle of incidence, reaching a maximum 
at 6 — 9' = \ , consistent with observations that indicate 
that a microtubule which is hindered in its growth will 
undergo a catastrophe at a rate that depends inversely 
on its growth speed [26] . 
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Figure 3: Probabilities for zippering, cross-overs and catas- 
trophes as deduced from the observations of [15] (combined 
data from MBD-DsRed and YFP-TUA6 labelling). Light grey 
shaded region: fraction of zippering events. Dark grey shaded 
region: fraction of induced catastrophes. White region: frac- 
tion of crossovers. Every data point is located at the centre 
of the corresponding bin, and the shaded regions have been 
extended to the boundaries using horizontal lines. The cor- 
responding lowest order Fourier coefficients of the interaction 
functions are: cq = 0.59, C2 = —0.36, zq = 0.24 (computed 
using numerical integration of the product of jsin(9| and a 
piecewise linear interpolation of the data). 



B. Continuum model 

Since there are many (« 10 2 - 10 3 ) microtubules, each 
of which can have multiple segments, in the cortical ar- 
ray of a typical interphase plant cell we treat the system 
using a coarse-grained description. In this approach, in- 
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IdDlC 1. 


Overview of all parameters and variables in natural dimensions 


Parameters 


Description 


Dimensions 


v+ 


growth speed 


[length] / [time] 


v~ 


shrinkage speed 


[length] / [time] 


r c 


catastrophe rate 


1/ [time] 


r r 


rescue rate 


1/ [time] 


In 


nucleation rate 


[length] ~ [time] " 1 


Pc(0) 


probability of induced catastrophe upon collision 1 




probability of zippering upon collision 


1 


Synthetic parameters 


^ V V >r 


growth parameter 


1/ [time] 


U = 1+ 2± 


speed ratio 


1 


c(0) =sin(|0|)P c (0) 


-> {c„} effective catastrophic collision probability 


1 


Z (0) = sin(|6»|)P z (6>) <— 


■> {in} effective zippering probability 


1 


Dependent variables 




microtubule length density 


[length] 1 [radian] 




average microtubule segment length 


[length] 


{m+(;,0),m-(Z,e),m 


% (1,9)} density of growing/shrinking/inactive segments 






with length / and direction 9 


[length] _ 3 [radian] _1 



stead of individual microtubules, we consider local den- 
sities of microtubule segments. This approximation is 
reasonable as long as the length scale of an individual 
microtubule segment is small compared to the linear di- 
mensions of the cell. From the outset we assume that the 
system is (and remains) spatially homogeneous, and we 
will eventually restrict ourselves to the steady-state solu- 
tions. In order to deal with the memory effect caused by 
the isotropic nucleation, followed by subsequent reorient- 
ing zippering events, we need to keep track of the segment 
number i, which starts at 1 for the segments connected 
to their nucleation site and increases by unity at each 
zippering event. Our fundamental variables are there- 
fore the areal number densities mf(l, 9, t) of segments in 
state a € {0, — , +} with segment number i having length 
I and orientation 9 (measured from an arbitrary axis) at 
time t. These densities obey a set of master equations 
that can symbolically be written as 

growth ~"T" ^rescue ^spont. cat. 
^induced cat. ^ zipper (1^) 
dtlTl^ {li, 9i,t) — ^shrinkage ^rescue ""T" ^spont. cat. 

~t~ ^induced cat. ~T" ^reactivation (-^) 
dt^i @i, t) zipper ^reactivation (1^) 

The flux terms $ event couple the equations for the grow- 
ing, shrinking and inactive segments and between differ- 
ent values of i. Equations (1) must be supplemented by 
a set of boundary conditions for the growing segments at 
I = 0. For the initial segment (i = 1) this reflects the 
isotropic nucleation of new microtubules, given by 

v+m+(h = 0,e,t) = ^, (2) 



where r n is nucleation rate. For the subsequent segments 
i > 1, this 'nucleation' of growing segments is the result 
of the zippering of segments with index i — I. Defining 
^zipper (Oi-i — * 0i,U-i,t) as the flux of i-segmcnts with 
angle 9i and length li zippering into angle 6i+i at time 
t (this will be made explicit in equation (13)), we obtain 
the boundary condition 

v+mf> 2 (l i = 0,6 i ,t) = 

J dli-i J d6i-i </?zipper (Oi-l — ► 6i, Z»_l,t) . (3) 

Generally, this leads to a qualitatively different bound- 
ary condition for every value of i. The model therefore 
consists of an infinite set of coupled equations, three for 
every value of i. However, in section II C we will show 
that in the steady state, this can be reduced to a finite set 
by summing over all segment indices i. In the following, 
we derive explicit expressions for each of the flux terms 

^ event • 

1. Growth and shrinkage terms: $ gr owth, ^shrinkage 

& growth m Equation (la) corresponds to the length in- 
crease of the growing segments. For segment growth in 
isolation, the length increase in a small time interval St 
is given by v + 5t, where v + is the growth velocity, and we 
have m + (l + v + St, 9,t + St) = m + (l, 9, t). By expanding 
the left hand term to first order in St, we find 

d t m+(l,9,t) = -v+dim+(l,9,t) = <S> growth (4) 

A similar derivation yields that 

d t m~(l,9,t) = v~dimj(l,6,t) = shrink (5) 
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where v is the shrinking velocity. 



2. Dynamic instability terms: ^rescue ,&s P ont. cat. 

^rescue and ® sp0 nt. cat. in equations (la) and (lb) cor- 
respond to the fluxes due to the spontaneous rescues 
and spontaneous catastrophe respectively and are simply 
given by 

^rescue =r r m^ (k, 9i,t) (6) 

$ spent, cat. =r c mf(li,6i,t) (7) 

where r r is the spontaneous rescue rate and r c is the 
spontaneous catastrophe rate. 

So far, we have described the first three terms of equa- 
tions (la) and (lb) (growth, shrinkage and dynamic in- 
stability terms). Together, these fully describe a system 
of non-interacting microtubules, in which also the bound- 
ary condition (3) vanishes due to the absence of zippering. 
In this special case we recover the well-known equations 
introduced by Dogterom and Leibler [22] (for i = 1). 



initions we can write the interaction terms as 

induced cat. = v+m+{k, 9 U t) J dO' c{9 % - 9')k(9' , t) 

(10) 

dipper = « + m+(ii,M) / d9' z(9 i -0')k{9',t) 

(11) 

The analogous term for crossovers is not used, because 
the occurrence of a crossover event has no effect on the 
growth of a microtubule. 

4- Reactivation term: $ reactivation 

The reactivation term ^reactivation corresponds to the 
flux of active microtubule segments, with segment num- 
ber i + 1, that, by shrinking to length zero, reactivate 
a previously inactive segment, effectively undoing a past 
zippering event and so create a new, shrinking, active 
segment with segment number i. The incoming flux of 
of such segments coming from a given direction 9i + i is 
given by v~ m~ +1 (li + i = 0, t). The reactivation flux 
is given by 



3. Interaction terms: ^induced cat. , zipper 

An interaction can occur when a growing active micro- 
tubule segment collides with another segment, irrespec- 
tive of the latter's state and length. This prompts the 
definition of the total length density k(9,t) of all micro- 
tubule segments in direction 9 at time t, given by: 

k{9,t)=J2 J d ^ k{mi{l h e,t)+mr{l i ,e,t)+m Q i {luO,t)) 

i 

(8) 

The density of collisions of a microtubule segment grow- 
ing in direction 9 with other segments in direction 9' is 
determined by the geometrical projection 

\sin(9-9')\k(9',t), (9) 

where the presence of the sin(# — 9') factor ensures the 
correct geometrical weighting reflecting the fact that par- 
allel segments do not collide. When a collision occurs, 
one of the three possible events, induced catastrophe 
(c) , zippering (z) or cross-over (x) occurs, with probabil- 
ities P c (9 - 8') , P z (0 - 9 1 ) and P x (9 - 6') respectively 
These probabilities can (and in-vivo do, see Figure 3) de- 
pend on the relative angle 9 — 9' between the incoming 
segment and the 'scatterer'. For convenience sake we 
absorb the geometrical factor |sin(0 — 9')\ into the prob- 
abilities, by defining / (0 - 9') = |sin(0 - 9')\ P f (9 - 8') 
for all events / €E {c, z, x} . The incoming flux of growing 
microtubule segments with given segment number, length 
and orientation is given by v + mf(l, 9, t). With these def- 



^rcactivation 

J d9 i+ i v~ml +l (k +1 = Q,8 t+ i,t)p nnziv (9iJ t \9 i+ i,t), 

(12) 

where the 'unzippcring' distribution Punzip(9i,h\9i+i,t) 
gives the probability that the shrinking microtubule reac- 
tivates an inactive segment with orientation 9i and length 
ij. This distribution will be determined below. 

A microtubule that has zippered will take a certain 
amount of time r to undergo a catastrophe and return 
to the zippering location, where r is a stochastic vari- 
able. The unzipppcring flux from direction Qi+i at time 
t consists of microtubules that had zippered at a range of 
times t — t and have now returned to the zippering loca- 
tion. This implicitly defines an originating time distribu- 
tion ^origin (t — 7"|0i+i,t) for the returning microtubules. 
Furthermore, because the evolution of a microtubule be- 
tween the zippering event and its return to the same lo- 
cation does not depend on the previous segments, the 
segment that is re-activated by a microtubule returning 
to the zippering position after a time r should be selected 
proportional to the 'forward' zippering flux at time t — r. 
The forward flux y z ip per (6i — * h, t) of microtubules 
with length li and angle 9i zippering into angle 9i + \ is 
defined in accordance with equation (11) as 

^zipper (9i — * Oi+1, k,t) = 

v + m+{k,8 l ,t)z(9 t - 9 l+1 )k(9 l+u t). (13) 

At each of the originating times t — r, the distribution 
of microtubules that zipper into the direction 9i + \ with 
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length li and orientation 9i is given by 

Vzippor (0 i ^ Qi+lji, t — t) 



Pzip(0i,k\0i+l,t-T) 



J dl'dO' <p* PP *(6'^0 i+1 ,l',t-T) 

(14) 

The probability distributions p rigin (t — r|#i+i, t) and 
Pzip(#i, h\0i+i, t — t) can be combined to determine the 
unzippcring distribution 



Punzip(0i,k\8i+l,t) = / dTPorigin (t - r\9 i+1 ,t) x 



^zipper {Oi — > ,k,t — 7") 

J dl'dO' <p zippet (0'^0 i+1 ,l',t-T) 



, (15) 



where we assume the system evolved from an initial con- 
dition at to = in which no microtubules were present. 
Clearly all the complicated history dependence of the sys- 
tem is hidden in the originating time distribution. How- 
ever, in the steady state situation we consider below, the 
time-dependence drops out and the details of this distri- 
bution become irrelevant. 



C. The steady state 

We now consider the steady state of the system of 
equations we have formulated. Setting the time deriva- 
tives to zero, the sum of equations (la) to (lc) yields 
^growth + ^shrinkage = 0, which together with equation 
(4) and (5) implies d lt [v + mf(li,9i)~v~mY(k,9i)] 
0. Because physically acceptable solutions should be 
bounded as U — > oo, we obtain the length flux balance 
equation 



v + mf(k,6i) = v m i (Z»,0»), 



(16) 



showing that the growing and shrinking segments have, 
up to a constant amplitude, the same oricntational and 
length distribution. This allows us to eliminate m~(li, 9i) 
from (la) to obtain 

di i mf{l i ,9i) = mf(k,9i)x 

g - [ d9' [c(|0 4 - 0'\) + z(\9i - 9'\)} k(9')\ , (17) 



(18) 



where the growth parameter 



characterizes the behaviour of the bare, non-interacting, 
system in which microtubules remain bounded in length 
for g < and become unbounded for g > 0. As the brack- 
eted factor on the right hand side of (17) does not depend 
on the segment length nor on the segment number, we 
immediately obtain that m^{li,9i) has an exponential 
length distribution 



where the average segment length / (9i) in the direction 
9i is given by 



g + / d9' (c(9 - 9') + z(9 - 9'))k{9'). (20) 



1 

W) 

The nucleation boundary conditions (2) and (3) are now 
transformed into independent nucleation equations that 
are expressed in terms of the amplitudes mf(9i) 



vrn{{9) 



2tt' 



(21) 



m+> 2 (6) = k{9) J d9' z{9' - 9)l(9')m+_ 1 (9'). (22) 

We now note that under these conditions equation (lc) 
is already satisfied, as can be explicitly checked by con- 
sidering ^reactivation (12) and using that in the steady 
state </J z ipp C r does not depend on time and the integral 
over ^origin is by definition equal to 1. This, in combi- 
nation with the results (16), (19) and (22), yields the 
identity with $ Z ip P or- We therefore need an independent 
argument to fix the densities of the inactive segments. To 
obtain this wc use the steady state rule that population 
size = nucleation rate x average lifetime. Consider a newly 
'born' growing segment, created cither by a nucleation or 
a zippcring event. Its average life time is by definition 
the average time until it shrinks back to zero length, i.e. 
the average return time. Clearly this time only depends 
on its orientation 9, the steady state microtubule length 
density k (9') and the dynamical instability parameters, 
but not on the segment number. We therefore denote it 
by t (9) . The steady state density of inactive segments 
with length li , orientation 9 and segment number i is then 
given by 



m° (h, 9) = / d9' ^ zippcr (9 -> 9', h) r (9') 



(23) 



where <^ z ip pe r is defined by equation (13), as inactive seg- 
ments are created by a zippcring event. Because the only 
length-dependent term on the right-hand side is m + (/j, 9), 
it follows that the length-dependence of the inactive seg- 
ment distributions is proportional to those of the active 
segments, i.e. 



(24) 



At the same time, the total integrated length density of 
segments, both active and inactive, with segment number 
i + 1 in the direction 9' is given by 



JVgt»i (9') = J dl i+1 x 
[<i (h + i,0') + m+ +1 (h +1 ,9') + mr +1 (l i+1 ,6')] 

d9" [ dk<p zippei (8" ^6',k)T(6'), (25) 



mf(k 



-u/m) 



(19) 



where the last equality follows from the fact that every 
segment with index i + 1 has been created by a zippcr- 
ing event of a segment with index i. We solve this for 
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t {9') and insert the result in (23), which, after expand- 
ing y'zipper (13), produces the following expression for 
m° (9): 



m° 2 (9)=mf(9) J d9' z(9 - 9')l(9')x 
[m° i+1 {9') + mi +1 (6') + 7 



H+l 



Jd9"z(9" -9')l(9")mf{9") 



(26) 



We use the nucleation equation (22) to replace the in- 
tegral in the denominator of the integrand on the right 
hand side of this expression. In addition, we define the 
quantity Qi(9) through 

m? (0) = Qi(0) [m+(9)+m-(9)} 



(27) 



l + —)Qi (0) m+ (9) 



uQ % [9) m+ (9) 



and equation (26) leads to the following recursion relation 
for Q, (9) 



Qi (9) = J d9' 



(9-9')k(9')l{9')(l + Q l+1 (9')). 



(28) 

We now argue that the ratio Qi (9) is in fact independent 
of the segment number. Using the fact that the growing, 
shrinking and inactive segments have an identical expo- 
nential profile, it follows from (27) that Qi(9) is equal to 
the ratio between inactive and active segments 



Qi(0) 



N?(9) 



umf (9) N+ (9) + Nr (9) 



(29) 



After a new microtubule segment has been created it 
will generally spend some time in an active state and 
some time time in an inactive state. The expected life- 
time t(9) can also be separated into the expected ac- 
tive and inactive lifetimes for any newly created segment: 
t (#) = ^active (6) + T inac tivc (9). These lifetimes are nec- 
essarily proportional to the total number of active and 
inactive segments, so that Qi(9) = T inactivc (9) /T act i VC (#)■ 
As we have argued before, these lifetimes do not depend 
on the segment number, and, hence, neither docs Qi{9). 
An alternative route to the same conclusion follows from 
expanding out the forward recursion in (28) to show that 
Qi (9) can for every i formally be written as the same in- 
finite series of multiple integrals involving z (9 — 9'), k (9) 
and I (9). We therefore write the self-consistency relation- 
ship 

Q(9) = [ d9' z(9- 9') k (9 r ) I (9') (1 + Q (9 1 )) . (30) 



steady state 
i J 

= ul{9f{l + Q{9))Y^<{0)- (31) 



D. Dimensional analysis 

In order to simplify our equations for further analysis 
and to identify the relevant control parameter we per- 
form a dimensional analysis. We therefore introduce a 
common length scale and rescale all lengths with respect 
to this length scale. For example, our primary variables 
m~l (9) have dimension [length] 3 [radian] 1 . Taking our 
cue from (21) and (31) we adopt the length scale 



1 v 

7T U 



2tt 



(32) 



where the additional factor of 7r _1 within the parentheses 
is added to suppress explicit factors involving it in the 
final equations. This definition allows us to define the 
dimcnsionlcss variables 



L(9) 
K(9) 
M+ (9) = ?m 
G = gl 



I (9) /l 

ixk (9) l 



(33a) 
(33b) 
(33c) 
(33d) 



In the absence of interactions, (20) shows that the av- 
erage length I of the microtubule is given by / = —1/g. 
This implies G = —Iq/1, meaning that, for G < 0, G 
can be interpreted as a measure for the non-interacting 
microtubule length. 

In addition, we adopt the dimcnsionlcss operator nota- 
tion 

1 f 27r 

F [h] {9) = - / d9' f{9- 9') h (9') , (33e) 
I" Jo 

where F e {C, Z}. We are now in a position to ex- 
press the equations in terms of the dimensionless quan- 
tities. Applying (33c) to the nucleation equations (21) 
and (22) yields expressions for the growing segment den- 
sities M^(9). Furthermore, the M^(9) for the different 
segment labels can be absorbed into a single microtubule 
plus end density (density of active segments), given by 



oo 



T{9) = uL{9)^Mt{9) 

i=i 



(34) 



The final closure of this set of equations is provided 
by the definition of the length density (8) applied to the 



Performing all substitutions, the final set of dimension- 
less equations reads 
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1 



L(0) 
K{0) 
Q(0) 



Segment length 
= —G + C [K] {0) 



Z[K] (0) 



Density 

--L(0)(1 + Q(0))T(0) 
Inactive-active ratio 
= Z[LK(1 + Q)\ {9) 
Plus end density 
T(0) = L (6) + L(0)K (0) Z [T] (0) 



(35a) 

(35b) 
(35c) 
(35d) 



with 



G = 



2v+v~ 



r n (v 



+ 



?, 




--) 




(- 


V+J 



--4 (35e) 



Looking at the resulting equations, we see that the 
segment length L is determined by the intrinsic growth 
dynamics (G) and the collisions leading to induced catas- 
trophes and zippering. The segment length density K 
is the product of the plus end density, the ratio of all 
segments to active segments (1 + Q) and the average seg- 
ment length. The ratio Q of inactive to active segments 
is modulated by the zippering operator, and the plus end 
density T consists of contributions from direct nucleation 
and zippered segments. We only consider parameter re- 
gions with physically realizable solutions that have have 
real and positive values for L, K, Q and T. 

Finally, we note that the interaction operators defined 
by (33e) are convolutions of the operand with the interac- 
tion functions c(0) and z(9). Both interaction functions 
are symmetric and 7r-periodic, and can therefore be writ- 
ten in terms of their Fourier coefficients as 



f(6) = 4 + E /2nCOs(2^) 
n=l 

d0 f(0)cos(2n0) 



fin — — 



(36) 



(37) 







Using the identity cos(0 — 0') = cos(6')cos(0 / ) + 
sin(0)sin(0') we find that the functions cos(2n0) and 
sin(2n#) are eigenfunctions of the operators C and Z, 
with the Fourier coefficients ci n and Z2n, respectively, as 
eigenvalues: 



F [cos(2n0)] = f 2n cos(2n(, 



(38) 



This convenient property will be exploited in later sec- 
tions. 



III. RESULTS 

A. Isotropic solution 

In the isotropic phase all angular dependence drops 
out. Because C[l] = cq and Z[l] = zq we are left with 



the set of equations 
1 

I " 
K = 

Q~- 
f = 



-G + (c + z ) K 

L(1 + Q)T 
z Q LR(l + Q) 
L + z LKT 



(39a) 

(39b) 
(39c) 
(39d) 



where the overbar denotes quantities evaluated in the 
isotropic phase. Solving for Q and T and inserting this 
into equation (39b) readily gives 



K 



(i - zoLRy ' 



(40) 



which can be combined with equation (39a) to yield the 
following relationship between G and the density 



K(c K-G) 2 = 1. 



(41) 



We see that the isotropic density is an increasing function 
of the microtubule dynamics parameter G and does not 
depend on the amount of zippering. This can be under- 
stood by the fact that zippering only serves to reorient 
the microtubules, which has no net effect in the isotropic 
state. In the absence of induced catastrophes (co = 0), 
the density K diverges as G | 0, consistent with the re- 
sult by Dogterom and Leibler [22]. In the presence of 
induced catastrophes a stationary isotropic solution ex- 
ists for all values of G, although this solution need not 
actually be stable. 



B. Bifurcation analysis 

We now search for a bifurcation point by considering 
the existence of steady-state solutions which are small 
perturbations away from the isotropic solution. These 
solutions are parametrized as follows 



L = L(1 + X) 
K = K (1 + k) 

Q = Q0- + X) 
T = T(1 + t) 



(42a) 
(42b) 
(42c) 
(42d) 



Inserting these expressions into (35), subtracting the 
isotropic solutions and expanding to first order in the 
perturbations gives 



A = -N (C[k] + Z[«]) 
K = X + t + z Nx 



1 



X 



— Z A + 

zq 

t = X + N(z q k + Z[t}) 



(43a) 
(43b) 

(43c) 

(43d) 



where N = LK. Note that in these equations, N has 
become the control parameter instead of G. Using (43b) 
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and exploiting the linearity of Z. we expand 

Z [k] = Z [t] + Z [A + k + z N X ] - Z [k] (44) 
= ^(t-X)-z k + z x-Z[k]. (45) 

Solving this for Z [k] and inserting the result into equa- 
tion (43a), combined with (43b), yields the relation 



(1 - z N)k = -2NC[k] 



(46) 



In the absence of induced catastrophes (C[k] = 0; only 
zippcring), a bifurcation can only occur if zqN = 1, which 
only happens for diverging density and G = (as seen 
from equations (40-41)), therefore in this case there is no 
bifurcation. In the generic case where induced catastro- 
phes are present, (46) can be satisfied only if k{9) is an 
eigenfunction of C. We know that the family of functions 
cos(2n#), n > 1, are eigenfunctions of C with eigenvalues 
c 2n , and therefore get a set of bifurcation values for N, 
one for each eigenvalue: N$ n = (— 2c2 n + zq)^ 1 . In addi- 
tion, we know that the isotropic solution must be stable 
as G — > —oo, because in this limit the microtubules have 
a vanishing length and do not interact. Therefore, the 
relevant bifurcation point is that for the lowest value of 
G, corresponding with the most negative eigenvalue of 
C (see also section HIE). Assuming that the induced 
catastrophe probability increases monotonically with the 
collision angle, £2 is always the most negative eigenvalue, 
so 



N* = 



1 



-2c 2 + z 



(47) 



We now derive the location of this bifurcation point in 
terms of the control parameter G. Denoting N = LR, 
equation (40) can be transformed to N(l — zqN) 2 = Z 3 , 
into which we can substitute GL = (c + Zq)N — 1 from 
equation (39a) and solve for G giving 

G 3 N(1 - z,qN) 2 = [(c + z )N - 1] 3 (48) 
Combining this with the result (47) yields 



G* = (-2c 2 ) 



1/3 



-2c 2 



(49) 



The implication is that the location of the bifurcation 
point as a function of the control parameter G is deter- 
mined entirely by the eigenvalues of the induced catastro- 
phe function c{9). Like the density in the isotropic phase, 
the location of the bifurcation point, this time perhaps 
more surprisingly, does not depend on the presence or 
amount of zippering. 



C. Segment length and mesh size 

An attractive interpretation of the microtubule length 
density K(9) is that it represents the density of 'obstacles' 



that are pointing in the direction 8 as seen by a micro- 
tubule growing in the perpendicular direction. From the 
obstacle density we can define a mesh size - the aver- 
age distance between obstacles. Taking into account the 
geometrical factor sin(0), we obtain 



i 

irl 



dO' \sin{6-9')\K(6') 



(50) 



In the case of the isotropic solution, this simplifies to 
£ = ttIq I (4iT) . Using this equality we can derive an 
expression for the average microtubule length A in the 
isotropic phase, expressed in units of the mesh size. The 
length of each segment is given by L and the number of 
segments per microtubule is given by (1 + Q), so using 
(40) we find 



A 



l L(l + Q) _ 4i? 3 / 2 



(51) 



Inserting this result into (41) provides the relationship 
between A and G 



G 



4 

7rA 



ttcqA 



1 



(52) 



As was the case for the density, we see that the micro- 
tubule length as a function of mesh size does not depend 
on the amount of zippering. However, it should be noted 
that the mesh size is defined through the average distance 
between single microtubules. In real systems, zippering 
would naturally lead to bundling, which in turn produces 
a system that has a larger mesh size between bundles (see 
also the discussion). 

Combining equations (51) and (47), the expression for 
A at the bifurcation point becomes 

A* = — — (53) 

7TC 2 

Assuming a monotonically increasing induced catastro- 
phe probability P c {9), we know that the minimum value 
for C2 is reached when every collision at an angle larger 
than 45° leads to a catastrophe. From (53), we see that 
this implies A* > 3/(2^2), meaning that for a bifurcation 
to occur, the microtubules need to be longer (sometimes 
much longer) than the mesh size, as is to be expected. 

Equation (52) can also provide an interpretation of the 
length scale Iq. In the absence of catastrophic collisions, 
we find 



Ak=o = -(-G- 3 ) = i (I 



(54) 



where I = — 1/g is the average length of the microtubules. 
Iq is therefore a measure of the microtubule length that 
is required to enable a significant number of interactions 
(A = 4/7T for I = Iq). If the free microtubule length I is 
(much) shorter than Zo, the system is dominated by the 
(isotropic) nucleations, keeping the system in an isotropic 
state. On the other hand, when / 3> Iq, the interactions 
dominate and, depending on the interaction functions, 
the system has the potential to align. 
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Figure 4: Bifurcation diagrams for the simplified interaction functions using three different induced catastrophe parameters. 
The figures on the left depict the probability P c (9) to induce a catastrophe upon collision, along with the corresponding values 
of Co. The centre and right columns depict the corresponding bifurcation diagrams as a function of G, expressed in terms of 
the total density -Ktotai and the 2D nematic order parameter S2 , respectively, where Jftotai = / K(8)dd. The isotropic solutions 
are by definition disordered, so S2 = 0, and their density is computed from (41). The bifurcation point is determined using 
(49), with £2 = —1/2. For each diagram, ordered solutions have been computed for zo = (black), zo = 1 (blue) and zq = 10 
(red). The solutions have been computed using the method discussed in appendix A. Solid lines indicate stable solutions and 
dashed lines indicate unstable solutions (see also section HIE). Note that the case of Co = 1 in the absence of zippering is a 
singular case where the stability cannot be determined, because non-isotropic solutions only exist for G = 0. This has been 
indicated by a dotted line. The SVdiagrams include the asymptotic limit point at G — with absolute ordering (at infinite 
density). The labels a, b and c indicate the parameter values of the results depicted in figure 5. The fact that the solutions for 
52 in the case Co = § do not reach the asymptotic point (G = 0, S2 ~ 1) is a consequence of the slowdown in convergence of 
the path-following method as G f 0. 



D. Ordered solutions for simplified interaction 
functions 



To find solutions beyond the immediate vicinity of 
the bifurcation point, we are hampered by the fact that 
these solutions are part of an infinite-dimensional solu- 
tion space. In appendix A it is shown that the solutions 
can be constrained to a finite-dimensional space by re- 
stricting the interaction functions c(9) and z(Q) to a finite 
number of Fourier modes. 

In this section, we will define a set of simplified interac- 
tion functions by restricting ourselves to Fourier modes 
up to and including cos(40). These modes provide us 
with just enough freedom for the model to exhibit rich 
behaviour. Using the fact that c(0) = z(0) = z(w/2) = 0, 
we find that £2=0 and that both £4 and £4 are deter- 
mined by the remaining parameters. Furthermore, we 



introduce an overall factor of a in both equations, allow- 
ing us to set £2 = —1/2, so that c(%/2) = a. We thus 
obtain a system that is fully specified by the parameters 
Co, zq and a. 



c(9) = a 
z(9) = a 



|-icos(20) + i(l-c o )cos(40) 



20 



-^(l-cos(40)) 



(55a) 
(55b) 



For a = 1, Co and z$ are the actual Fourier coeffi- 
cients of the interaction functions. Demanding that 
Pc{9) = c(8)/sm(8) is monotonically increasing on the 
interval [0,7r/2] leads to the constraint 



3 „ 9 
4-°-8 



(56) 
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and zo is a positive real number. Of course, the total 
probability of zippering and catastrophe induction may 
not exceed 1, placing an upper bound on a. In the ab- 
sence of zippering, we have a < 1. 

It should be noted that the value of a has no quali- 
tative effect on the results. This can be understood by 
realizing that the set of equations (35) is invariant under 
the substitutions 



a) 



20 



C -> aC 
Z^aZ 

G -► a 1/3 G 



L 
K 
T 



-1/3 L 
-2/3 K 
-1/3 T 



where the first substitution reflects the presence of a in 
the definitions (55). The existence of this scaling relation 
implies that all functional dependencies between any of 
these parameters and variables remain unchanged when 
subjected to the inverse scaling. Explicitly, the relevant 
parameters become C/a, Z/a and a~ x ^G and the vari- 
ables a 1/3 L, o?^K and a 1/3 T. With this in mind, we 
have used the arbitrary choice a = 1 for our numerical 
calculations, indicating the appropriate scaling on the 
axes of figures 4 and 5. 

Equation (49) indicates that, for the simplified inter- 
action functions, the bifurcation point is located in the 
range 



- < G* < - 
4 " " 8 



(57) 



and from (41) and (53) we find that K* = oT 2 ! 3 and 
A* = 4/(o!7r). We have used the numerical procedure 
described in appendix A to determine the ordered solu- 
tions of (35), starting from the bifurcation point. This 
has been done for nine different parameter values. For 
the values of Co we used the extreme values 3/4 and 9/8, 
as well as 1, the latter corresponding to G* = 0. For each 
of these three cases, we have varied the zippering param- 
eter Zo, choosing values of 0, 1 and 10. Figure 4 shows 
the results, depicting both the total density of the sys- 
tem and the degree of ordering as a function of G. The 
degree of ordering is measured by the order parameter 
S2, defined as 



S 2 [K(0)} 



Jo n d6K(6) 



(58) 



which is the standard 2D ncmatic liquid crystal order 
parameter (yielding for a completely disordered system 
and 1 for a fully oriented system). 



E. Stability of solutions 

The bifurcation constraint (46) indicates that the space 
of bifurcating functions n n is spanned by the functions 
cos(2n#) and sin(2n6>) for a given value of n > 1. These 
solutions are therefore symmetric with respect to an ar- 
bitrary axis that we choose to place at 8 = 0. Even after 



G/oc 1/3 =-0.2 
z =0 



x 2 ^K 



a 



1/3 T 



-n/1 n/1 -n/1 n/1 



b) 



20 



G/oc 1/3 =-0.15 
z =0 



a 2 ' 3 K 





-n/1 n/1 -n/1 n/1 



c) 



20 




G/a 1/3 =-0.15 
z =l 



Figure 5: Three stable ordered solutions that correspond to 
the points labelled (a) , (b) and (c) in figure 4. The (unstable) 
isotropic solutions for the same parameter values are indicated 
with a dashed line. The parameter values for (a) and (b) differ 
only in the value of G, whereas the parameter values for (b) 
and (c) differ only in the value of zo. All results have been 
calculated using the method described in appendix A. 



the restriction to this symmetry axis there are still two 
solution branches emanating from the bifurcation point, 
differing in the sign of the coefficient of the perturbation. 
These branches correspond to solutions peaked around 
8 = and 8 = 7r/(2n), respectively, that arc identical ex- 
cept for this rotation. The symmetry of these solutions 
indicates that the bifurcation is of the pitchfork type. In 
figure 5 we plot the solutions with a maximum at 8 = 0. 

The presence of a pitchfork bifurcation implies a loss 
of stability of the originating branch [27]. In our sys- 
tem, we know that the isotropic solution must be stable 
in the limit G — > —00. Therefore, the local stability 
of the isotropic solution is lost at the first bifurcation 
point (for the lowest value of G), corresponding to the 
eigenfunction cos(2#). Because this eigenfunction is or- 
thogonal to the eigenfunctions related to the subsequent 
bifurcation points (005(2716*), n > 2), the stability of the 
unstable mode will not be regained at any point along the 
isotropic solution and the isotropic solution itself remains 
unstable for all higher values of G. This also means that 
the solution branches originating at further pitchfork bi- 
furcations will be unstable near the isotropic solution. In 
this paper we restrict ourselves to the analysis of the first 
bifurcation point and the corresponding ordered solution 
branch. Because the solutions on this branch already 
have the lowest symmetry permitted by the interaction 
functions, there arc no further bifurcation points along 
this branch. 
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Generically [27, chap. IV], the branches of the initial 
pitchfork bifurcation are stable for a supercritical bifurca- 
tion (branches bending towards higher values of G) and 
unstable for a subcritical bifurcation (branches bending 
towards lower values of G). In addition, turning points 
in the bifurcating branches generally correspond to an ex- 
change of stability [28, page 22]. This analysis allows us 
to assign stability indicators to the bifurcation diagrams 
in figure 4, even in the absence of a detailed study of the 
time-dependent equations (1). 



IV. DISCUSSION AND CONCLUSIONS 

Based on biological observations, we have constructed 
a model for the orientational alignment of cortical micro- 
tubules. The model has a number of striking features. 
First of all, for a given set of induced catastrophe and 
zippering probabilities (P c {0) and P z (9)), it allows us 
to identify a single dimensionless control parameter G, 
which is fully determined by the nucleation rate and in- 
trinsic dynamics of individual microtubules. This result 
by itself may turn out to be very useful in comparing 
different in vivo systems or the same system under differ- 
ent conditions or in different developmental stages. For 
increasing values of G, the isotropic stationary solutions 
to the model show an increase both in density and in 
abundance of interactions, as measured by the ratio of 
microtubule length to the mesh-size. Secondly, the bifur- 
cation point, i.e. the critical value of G* of the control 
parameter at which the system develops ordered station- 
ary solutions from the isotropic state, is determined solely 
by the probability of collisions between microtubules that 
lead to an induced catastrophe. 

Indeed, from the numerical solutions of the minimal 
model introduced in Section HID it appears, perhaps 
surprisingly, that the co-alignment of microtubules due 
to zippering events, if anything, diminishes the degree 
of order. These results identify the "weeding out" of 
misaligned microtubules — by marking them for early 
removal by the induced switch to the shrinking state — 
as the driving force for the ordering process. Finally, 
in spite of not being able to directly assess the stability 
of the solutions in the time-domain, we have provided 
arguments that stable ordered solutions are possible for 
the regime G < 0, i.e. where the length of individual 
microtubules is intrinsically bounded. 

For G > 0, individual microtubules have the tendency 
to grow unbounded, unless they are kept in check by 
catastrophic collisions. Although (locally) stable solu- 
tions may exist for values of G that are not too large, for 
every G > there exists a class of aligned 'runaway' so- 
lutions with diverging densities. The computed ordered 
solutions, regardless of their stability, converge to a point 
with G = 0, for which the microtubules are perfectly 
aligned (S2 = 1) and the system is infinitely dense. The 
existence of this point can be understood by the fact that 
the alignment also serves to decrease the number of colli- 



sions, and in the limit of a perfectly aligned system, the 
(relative) number of collisions vanishes. 

How realistic is the model presented? To answer 
this question we need to consider several known factors 
that have not been included. First of all, microtubules 
typically can de-attach from their nucleation sites and 
then perform so called treadmilling motion, whereby the 
minus-end shrinks at a more or less steady pace, which 
is small compared to both the growth- and the shrinking 
speed of the more active plus end. In the case that no zip- 
pering occurs at all it is relatively easy to show that the 
effect of treadmilling simply leads to a rcnormalization of 
the parameter G and the interaction functions c(9) and 
z(9), but leaving the qualitative behaviour of the model 
identical to the one discussed here. When zippering does 
occur, one expects the treadmilling to enhance the degree 
of ordering in an ordered state, as over time it "eats-up" 
the, by definition less ordered, initial segments of each 
microtubule. This effect is also consistent with the ob- 
servation in figure 5c that in the case with zippering the 
active tips are on average more strongly aligned than the 
average segment. In fact, given that the comparison be- 
tween figures 5b and 5c also shows that, all else being 
equal, zippering sharpens the orientational distribution 
of the active tips as compared to the case with no zipper- 
ing, it is conceivable that the combination of zippering 
and treadmilling could lead to more strongly ordered sys- 
tems for the same value of the control parameter. 

Next it is known that in vivo severing proteins, such 
as katanin are active in, and crucial to, the formation of 
the cortical array [29] . Although in principle the effect of 
severing proteins could be included in the model, it would 
present formidable problems in the analysis as well as 
introduce additional parameters into the model for which 
precise data is lacking. 

Another effect that has not been taken into account 
explicitly is microtubule bundling. Whenever a micro- 
tubule zippers alongside another segment, they form a 
parallel bundle [30]. However, the coarse-grained nature 
of our model precludes the formation of bundles and only 
allows for alignment of the segments. This means that 
a microtubule that is growing in a different direction en- 
counters each microtubule separately rather than as a 
single bundle. It is to be expected that the catastrophe 
and zippering rates stemming from N individual colli- 
sions will be higher than those from a single collision 
with a bundle of N microtubules. Hence, in realistic 
systems the event rate is likely to be lower than that pre- 
dicted by the model, or, in other words, corresponds to 
a lower value of the scaling parameter a. Bundles are 
also thought to be more than just adjacently aligned mi- 
crotubules, because they may be stabilized through as- 
sociation with bundling proteins that could potentially 
decrease the catastrophe rate of individual microtubules 
within a bundle (sec [31]). This is a non-trivial effect that 
should be considered separately and is likely to lead to 
an increased tendency to form an ordered structure. 

Finally, our model implicitly assumes that there is an 
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infinite supply of free tubulin dimers available for incor- 
poration into microtubules. Although there is no defi- 
nite experimental evidence for this, it is reasonable to 
assume that in vivo there is a limit to the size of the 
free tubulin pool. Such a finite tubulin pool would have 
marked consequences for the behaviour of the model, 
because the growth speed, and possibly also the nuclc- 
ation rate, are dependent on the amount of free tubu- 
lin, or equivalently the total density of microtubules ktot- 
To a first approximation the growth speed is given by 

v + (k to t) = v + (k tot = 0) (l - jr^J where k max the max- 
imally attainable density when all tubulin is incorporated 
into microtubules. This allows for stable states to de- 
velop even when G(k to t = 0) > 0, because under this 
pool-size constraint the length of individual microtubules 
will always remain bounded, and the system will settle 
into a steady state with G(ktot) < G(k to t = 0). This 
behaviour could provide a biologically motivated mech- 
anism by which a solution with a particular density is 
selected. 

To see whether our model, in spite of its approximate 
nature, makes sense in the light of the available data we 
first use the collision event probabilities obtained by Dixit 
and Cyr [15] (see figure 3) to obtain an estimate for the 
bifurcation value of the control parameter of G* = —0.15 
for the case of Tobacco BY-2 cells. An ordered phase 
of cortical microtubules should therefore be possible pro- 
vided G > G*. Given the available data on the micro- 
tubule instability parameters in this same system taken 
from Dhonukshc et al. [10] and Vos et al. [32] we would 
predict using the definition (35e) that this requires the 
nucleation rate of new microtubules to be larger than 0.05 
min _1 /xm -2 (Dhonukse) and 0.01 min" 1 /! m~ 2 (Vos) re- 
spectively. Both these estimates for a lower bound on 
the nucleation rate are reasonable as they imply the nu- 
cleation of order 10 3 microtubules in the whole cortex 
over the course of the build-up towards full transverse 
order, comparable to the number that is observed. 

Finally, we should point out that our model so far 
only addresses the question of what causes cortical mi- 
crotubules to align with respect to each other. Given 
that in growing plant cells the cortical array is invariably 
oriented transverse to the growth direction, the question 
of what determines the direction of the alignment axis 
with respect to the cell axes is as, if not more, important 
from a biological perspective. We hope to address this 
question, as well as the influence of some of the as yet 
neglected factors mentioned above, in future work. 
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Appendix A: NUMERICAL EVALUATION OF 
THE ORDERED SOLUTIONS 

The solutions to the set of equations (35) lie in an 
infinite-dimensional solution space. This creates signifi- 
cant hurdles for the numerical search for solutions. In 
this section, we will see that it is possible to restrict 
the solutions to a finite-dimensional space by imposing 
constraints on the interaction operators C and Z. In 
addition, we present a method to follow the branch of or- 
dered solution in this finite-dimensional space, starting 
from the bifurcation point (49) . 

We start by reformulating the set of equations (35) by 
replacing L(6) and T(0) through the definitions 



S(6) 



1 



L(0)> 



U{0) 



1 



( T{9) 
K{0) U(0) 



1 



(Al) 



Following these substitutions, the interaction operators 
are all applied at the outermost level of the equations, 
enabling us to make use of their properties in Fourier 
space. Explicitly, we obtain 



and 



s(0) 
Q{0) 

U{0) 



K{0) 



-G + C [K] {9) + Z [K] {0) 
Z[K(1 + Q)/S] (0) 
Z[(l + KU)/S] (0) 



1 + Q(0) 
S2(0)-U(0)(1 + Q(0)Y 



(A2) 
(A3) 
(A4) 



(A5) 



Denoting the Fourier components of S(8), Q{0) and U (0), 
by s n , q n and u n , respectively, the interacting micro- 
tubule equations reduce to a (potentially infinite) set of 
scalar integral equations: 



S2n — —28n.oG 



Cln + Z2r, 



2ji 



d9 cos{2n9)K(9) 



Z2n 



hn = — / d0 



2jt 



Z2n 



U 2n = / d0 



2jt 



cos(2n0)K(0) (1 + Q(0)) 
W) 

cos(2n0)(l + K{0)U{0)) 
W) ' 



(A6a) 
(A6b) 



(A6c) 



From the structure of these equations, wc immediately 
see that we can greatly reduce the dimensionality of the 
problem by setting a number of Fourier coefficients Z2n 
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and c 2n to zero. In other words, by restricting our space 
of interaction functions c(9) and z(0), the problem can 
be reduced to a finite number of scalar equations. 

Applied to the simplified interaction functions intro- 
duced in section HID, we know that the sets of station- 
ary solutions form lines in the 8-dimensional phase space 
spanned by the variables {sq, §2, S4, <?o, 54, uq, U4} and the 
parameter G. At least two of such solution lines exist, 
one corresponding to the isotropic solution and the other 
to the ordered solution, and these lines intersect at the 
bifurcation point. 

Within this 8-dimensional space, we have used a nu- 
merical path-following method similar to the one de- 
scribed in [33, 34] that follows the ordered solution 
branch by searching for a local minimum in the root mean 
error of the constituent equations (A6). The search for or- 
dered solutions is initiated at the bifurcation point, with 



coordinates 

_ z - 2c 2 = zo_ _ zq 

(-2c 2 )2/3' v ■ 2c 2 ' (-202)1/3' 

(A7) 

so that in the case of our simplified interaction model 

{G, s , s 2 , S4, Qo, 94, "0, "4)0 = 

{c - 1, 2(z + 1), 0, 0, 22 , 0, 2z , 0}. (A8) 

The initial instability affects only the cos(26*) mode. This 
mode only appears in the equation for §2 and the remain- 
ing parameters are affected only by higher order correc- 
tions. For this reason we choose the initial direction of 
the path to be the unit vector in the ^-direction and the 
path is traced from there. 
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